The Data.

City-level information pulled from the following sources:

LA Almanac (race, age, and income) http://www.laalmanac.com/population/po38.php http://www.laalmanac.com/employment/em12.php

LA County vaccine percent by city http://publichealth.lacounty.gov/media/coronavirus/vaccine/vaccine-dashboard.htm

LA County age adjusted incidence rates http://dashboard.publichealth.lacounty.gov/covid19_surveillance_dashboard/

Data are incidence rates (proportions), vaccination rates (proportions), and characteristics by city.

Hypotheses:

All hypothesis are on the city-level and care must be taken when interpreting the results as they describe city-level variation and not individual level variation or co-variation.

Investigation of the December 2020-January 2021 peak for COVID 19 captured by a single day for incidence rates: 2021-01-07:

How do race/ethnicity, age, and income describe variation in incidence rates?

Investigation of the August 2021 peak for COVID 19 captured by a single day for incidence rates: 2021-08-18:

How do August vaccination proportions, race/ethnicity, age, income, and incidence rates for the Dec.-Jan. peak describe variation in incidence rates during the August peak for COVID 19?

Investigation of August 2021 vaccination rates for COVID 19 captured by a single day for vaccination rates: 2021-08-18:

How do race/ethnicity, age, income, and incidence rates for the Dec.-Jan. peak describe variation in August vaccination rates?

Variable Distributions and Coding

Cities

d <- readRDS("data/analysis_data.rds") %>% 
  arrange(Date.x) %>% 
  filter(Date.x == "2021-08-18")

city.populations <- data.frame(d$city, d$population.x)
names(city.populations) <- c("City", "Population_Size")
city.populations <- city.populations[order(city.populations$Population_Size, decreasing=T), ]
kable(city.populations, format = "simple", align = 'c', caption = "Cities and Population Size")

Cities and Population Size
	City	Population_Size
121	Los Angeles	4043337
95	Santa Clarita	220424
44	Glendale	206493
62	Lancaster	161570
78	Palmdale	158968
82	Pomona	157869
108	Torrance	149268
35	East Los Angeles	125269
39	El Monte	117269
33	Downey	114263
52	Inglewood	113582
114	West Covina	108235
77	Norwalk	107622
18	Burbank	107180
25	Compton	99904
101	South Gate	98155
20	Carson	93846
97	Santa Monica	92446
48	Hawthorne	91301
119	Whittier	91216
4	Alhambra	86724
61	Lakewood	80362
15	Bellflower	77735
12	Baldwin Park	76769
68	Lynwood	72047
85	Redondo Beach	68697
11	Azusa	65963
26	Covina	65851
6	Arcadia	65735
42	Florence-Graham	64705
74	Montebello	64375
81	Pico Rivera	64284
75	Monterey Park	62262
43	Gardena	61310
51	Huntington Park	59484
104	South Whittier	59222
32	Diamond Bar	57515
80	Paramount	56023
46	Hacienda Heights	55926
88	Rosemead	55350
45	Glendora	52764
89	Rowland Heights	51022
22	Cerritos	50067
56	La Mirada	49599
5	Altadena	43620
14	Bell Gardens	43071
84	Rancho Palos Verdes	42747
73	Monrovia	42681
8	Westmont	42442
92	San Gabriel	40954
57	La Puente	40697
29	Culver City	39865
115	West Hollywood	36951
23	Claremont	36484
107	Temple City	36455
13	Bell	36332
70	Manhattan Beach	35999
58	La Verne	35322
120	Willowbrook	34913
16	Beverly Hills	34520
90	San Dimas	34516
63	Lawndale	33614
111	Walnut	30532
72	Maywood	28049
21	Castaic	27191
34	Duarte	26444
102	South Pasadena	26053
91	San Fernando	24612
28	Cudahy	24347
19	Calabasas	24323
76	East San Gabriel	24036
110	Valinda	23371
100	South El Monte	22680
64	Lennox	22542
113	West Carson	22086
105	Stevenson Ranch	20966
2	Agoura Hills	20883
67	Lomita	20729
54	La Crescenta-Montrose	19801
49	Hermosa Beach	19670
96	Santa Fe Springs	18364
7	Artesia	16795
40	El Segundo	16786
112	Walnut Park	16143
37	East Rancho Dominguez	15308
47	Hawaiian Gardens	14676
79	Palos Verdes Estates	13522
93	San Marino	13277
27	Charter Oak	13144
24	Commerce	13069
60	Lake Los Angeles	12994
69	Malibu	12961
83	Quartz Hill	12906
99	Signal Hill	11797
98	Sierra Madre	10989
116	West Puente Valley	9835
71	Marina del Rey	9411
103	South San Gabriel	8848
118	Westlake Village	8360
87	Rolling Hills Estates	8113
1	Acton	7971
59	Ladera Heights	7071
10	Avocado Heights	6775
36	East Pasadena	6403
106	Sun Village	6036
55	La Habra Heights	5455
38	East Whittier	5306
30	Del Aire	4393
3	Agua Dulce	4158
66	Littlerock	4021
9	Avalon	3869
109	Val Verde	3309
31	Desert View Highlands	2493
94	San Pasqual	2035
86	Rolling Hills	1940
50	Hidden Hills	1890
65	Leona Valley	1751
41	Elizabeth Lake	1661
53	Irwindale	1459
117	West Rancho Dominguez	1359
17	Bradbury	1069

City <- d$city

Incidence Rates

January Peak: 2021-01-07

Jan.IR <- d$adj_case_14day_rate.y
Jan.IR.cat <- relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
                      ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(Jan.IR, breaks=seq(from=0, to=round(max(Jan.IR),-3), by=200), col="red", main="January Incidence Rates for Cities", xlab="Incidence Rates")
abline(v=quantile(Jan.IR, probs=c(.25, .5, .75)), lwd=2, lty=2)

August Peak: 2021-08-18

Aug.IR <- d$adj_case_14day_rate.x

hist(Aug.IR, breaks=seq(from=0, to=round(max(Aug.IR),-2), by=100), col="red", main="August Incidence Rates for Cities", xlab="Incidence Rates")

Percent Vaccinated in Each City

Aug.Vax.percent <- d$dose1_all_c_prcent
Aug.Vax.cat <- relevel(as.factor(ifelse(Aug.Vax.percent > quantile(Aug.Vax.percent, probs=.75), "High",
                      ifelse(Aug.Vax.percent < quantile(Aug.Vax.percent, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(Aug.Vax.percent, breaks=seq(from=0, to=100, by=5), col="blue", main="Percent Vaccination Rates for Cities", xlab="Percent Vaccinated")
abline(v=quantile(Aug.Vax.percent, probs=c(.25, .5, .75)), lwd=2, lty=2)

Race/Ethnicity

Proportion of Hispanic Individuals

CityHispanic <- d$Hispanic
CityHispanic.cat <- relevel(as.factor(ifelse(d$Hispanic > quantile(d$Hispanic, probs=.75), "High",
                      ifelse(d$Hispanic < quantile(d$Hispanic, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(CityHispanic, breaks=seq(from=0, to=100, by=5), main="Proportion of Hispanic Individuals", xlab="Proportion")
abline(v=quantile(CityHispanic, probs=c(.25, .5, .75)), lwd=2, lty=2)

Proportion of Black Individuals

CityBlack <- d$Black
CityBlack.cat <- relevel(as.factor(ifelse(d$Black > 10, "High","Low")), "Low")


hist(CityBlack, breaks=seq(from=0, to=100, by=5), main="Proportion of Black Individuals", xlab="Proportion")
abline(v=10, lwd=2, lty=2)
text(10, 60, "Cutoff at 10%", pos=4)

Proportion of Asian Individuals

CityAsian <- d$Asian
CityAsian.cat <- relevel(as.factor(ifelse(d$Asian > quantile(d$Asian, probs=.75), "High",
                      ifelse(d$Asian < quantile(d$Asian, probs=.5), "Low",
                             "Med"))), ref="Low")

hist(CityAsian, breaks=seq(from=0, to=100, by=5), main="Proportion of Asian Individuals", xlab="Proportion")
abline(v=quantile(CityAsian, probs=c(.25, .5, .75)), lwd=2, lty=2)

Age

Age data is captured by the number/proportion of individuals within each city in a age category. Age categories include: “Age.Under.15”,“Age.15.17”,“Age.18.24”,“Age.25.34”,“Age.35.54”,“Age.55.64”,“Age.65.”

Three groups are created for the analysis:

Young: combined age categories for “Age.Under.15”,“Age.15.17”,“Age.18.24”
Middle: age category for “Age.35.54”
Old: combined age categories for “Age.55.64”,“Age.65.”

Age <- d[,c("Age.Under.15","Age.15.17","Age.18.24","Age.25.34","Age.35.54","Age.55.64","Age.65.")]
Age <- t(apply(Age, 1, FUN=function(v) { v/sum(v) }))

CityAgeYoung <- apply(Age[,c("Age.Under.15","Age.15.17","Age.18.24")], 1, sum)
CityAgeOld <- apply(Age[,c("Age.55.64","Age.65.")], 1, sum)
CityAgeMid <- Age[,"Age.35.54"]
CityAge <- cbind(CityAgeYoung, CityAgeOld)

CityAgeOld.cat <- relevel(as.factor(ifelse(CityAgeOld > quantile(CityAgeOld, probs=.75), "High","notHigh")), "notHigh")
CityAgeYoung.cat <- relevel(as.factor(ifelse(CityAgeYoung > quantile(CityAgeYoung, probs=.75), "High","notHigh")), "notHigh")

ftable(CityAgeOld.cat, CityAgeYoung.cat)

##                CityAgeYoung.cat notHigh High
## CityAgeOld.cat                              
## notHigh                              61   30
## High                                 30    0

Household Income

HouseholdIncome <- d$Households
HouseholdIncome.cat <- relevel(as.factor(ifelse(HouseholdIncome > quantile(HouseholdIncome, probs=.75), "High",
                      ifelse(HouseholdIncome < quantile(HouseholdIncome, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(HouseholdIncome, breaks=seq(from=0, to=round(max(HouseholdIncome),-2), by=10000), main="Household Income Within Each City", xlab="Income", col="green")
abline(v=quantile(HouseholdIncome, probs=c(.25, .5, .75)), lwd=2, lty=2)

Variables for Analysis

Distributions

CityAge.m <- cbind(CityAgeOld.cat, CityAgeYoung.cat)
CityRaceEthnicty.m <- cbind(CityHispanic.cat, CityBlack.cat, CityAsian.cat)
d.analysis <- data.frame(Aug.IR, Jan.IR.cat, Aug.Vax.cat, Aug.Vax.percent, CityRaceEthnicty.m, CityAge.m, HouseholdIncome.cat)

summarytools::view(dfSummary(d.analysis, style = 'grid',
                               max.distinct.values = 10, plain.ascii = FALSE, valid.col = FALSE, headings = FALSE), method = "render")

Variable

Stats / Values

Freqs (% of Valid)

Graph

Missing

Aug.IR [numeric]

Mean (sd) : 436.8 (178.7) min < med < max: 39 < 447 < 1083 IQR (CV) : 186 (0.4)

111 distinct values

0 (0.0%)

Jan.IR.cat [factor]

1. Med 2. High 3. Low

61	(	50.4%	)
30	(	24.8%	)
30	(	24.8%	)

0 (0.0%)

Aug.Vax.cat [factor]

1. Med 2. High 3. Low

62	(	51.2%	)
29	(	24.0%	)
30	(	24.8%	)

0 (0.0%)

Aug.Vax.percent [numeric]

Mean (sd) : 69.8 (11.6) min < med < max: 20.9 < 70.6 < 91.8 IQR (CV) : 13.9 (0.2)

100 distinct values

0 (0.0%)

CityHispanic.cat [integer]

Mean (sd) : 1.7 (0.8) min < med < max: 1 < 1 < 3 IQR (CV) : 1 (0.5)

1	:	61	(	50.4%	)
2	:	30	(	24.8%	)
3	:	30	(	24.8%	)

0 (0.0%)

CityBlack.cat [integer]

Min : 1 Mean : 1.1 Max : 2

1	:	104	(	86.0%	)
2	:	17	(	14.0%	)

0 (0.0%)

CityAsian.cat [integer]

Mean (sd) : 1.8 (0.8) min < med < max: 1 < 2 < 3 IQR (CV) : 2 (0.5)

1	:	60	(	49.6%	)
2	:	30	(	24.8%	)
3	:	31	(	25.6%	)

0 (0.0%)

CityAgeOld.cat [integer]

Min : 1 Mean : 1.2 Max : 2

1	:	91	(	75.2%	)
2	:	30	(	24.8%	)

0 (0.0%)

CityAgeYoung.cat [integer]

Min : 1 Mean : 1.2 Max : 2

1	:	91	(	75.2%	)
2	:	30	(	24.8%	)

0 (0.0%)

HouseholdIncome.cat [factor]

1. Med 2. High 3. Low

61	(	50.4%	)
30	(	24.8%	)
30	(	24.8%	)

0 (0.0%)

Generated by summarytools 0.9.9 (R version 4.1.0)
2021-11-28

Correlation

# need to make sure low is lowest when setting as.numeric (to make cor more clear)





d.num <- data.frame(Aug.IR, Jan.IR.cat, Aug.Vax.cat,
                    CityHispanic.cat, as.numeric(CityBlack.cat), as.numeric(CityAsian.cat), 
                    as.numeric(CityAgeOld.cat), as.numeric(CityAgeYoung.cat),
                    HouseholdIncome.cat) %>% 
  mutate(across(c(Jan.IR.cat, Aug.Vax.cat, CityHispanic.cat, HouseholdIncome.cat), 
                ~ as.numeric(relevel(.x, ref = "Low"))))

cormat <- cor(d.num, use="complete.obs", method="pearson")
corrplot(cormat, type="upper", order="original",
          col=brewer.pal(n=8, name="RdYlBu"),
          title = "",
          addCoef.col = "black",
          tl.cex=.5, number.cex=.5)

Pairwise Trend Plots

January IR vs. Proportion of Hispanic Individuals

CityHispanic.cont <- as.numeric(relevel(CityHispanic.cat, "Low"))
coef <- summary(lm(Jan.IR~CityHispanic.cont))$coef
ggplot(d.analysis, aes(x=CityHispanic.cont, y=Jan.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=3000, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=2500, label=paste("p-value", "=", round(coef[2,4],25))) +
  scale_x_continuous(name ="Proportion of Hispanic Individuals", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("Jan Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

January IR vs. Household Income

coef <- summary(lm(Jan.IR~as.numeric(relevel(HouseholdIncome.cat, "Low"))))$coef
ggplot(d.analysis, aes(x=as.numeric(relevel(HouseholdIncome.cat, "Low")), y=Jan.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=3000, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=2500, label=paste("p-value", "=", round(coef[2,4],23))) +
  scale_x_continuous(name ="Household Income", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("January Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

August IR vs. January IR

coef <- summary(lm(Aug.IR~Jan.IR))$coef
ggplot(d.analysis, aes(x=Jan.IR, y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=3000, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=3000, y=800, label=paste("p-value", "=", round(coef[2,4],6)))+
  xlab("January Incidence Rates")+ylab("August Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

August IR vs. Proportion of Hispanic Individuals

CityHispanic.cont <- as.numeric(relevel(CityHispanic.cat, "Low"))
coef <- summary(lm(Aug.IR~CityHispanic.cont))$coef
ggplot(d.analysis, aes(x=CityHispanic.cont, y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=800, label=paste("p-value", "=", round(coef[2,4],2))) +
  scale_x_continuous(name ="Proportion of Hispanic Individuals", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("August Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

August IR vs. Household Income

coef <- summary(lm(Aug.IR~as.numeric(relevel(HouseholdIncome.cat, "Low"))))$coef
ggplot(d.analysis, aes(x=as.numeric(relevel(HouseholdIncome.cat, "Low")), y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=800, label=paste("p-value", "=", round(coef[2,4],2))) +
  scale_x_continuous(name ="Household Income", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("August Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

August IR vs. August Percent Vaccinated

coef <- summary(lm(Aug.IR~Aug.Vax.percent))$coef
ggplot(d.analysis, aes(x=Aug.Vax.percent, y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=80, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=80, y=800, label=paste("p-value", "=", round(coef[2,4],2)))+
  xlab("August Vaccination Percent")+ylab("August Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

January Incident Rates:

Regression Model

reg.Jan.noVax <- lm(Jan.IR ~ HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)

Regression summary

kable(summary(reg.Jan.noVax)$coef, digits=3, format = "simple", align = 'c', caption = "January Regression")

January Regression
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	1084.576	432.201	2.509	0.014
HouseholdIncome.catHigh	-982.530	171.650	-5.724	0.000
HouseholdIncome.catLow	784.904	164.034	4.785	0.000
CityRaceEthnicty.mCityHispanic.cat	96.811	82.005	1.181	0.240
CityRaceEthnicty.mCityBlack.cat	64.883	166.059	0.391	0.697
CityRaceEthnicty.mCityAsian.cat	68.113	73.316	0.929	0.355
CityAge.mCityAgeOld.cat	-106.842	147.045	-0.727	0.469
CityAge.mCityAgeYoung.cat	446.460	165.292	2.701	0.008

#stargazer(reg.Jan.noVax, type = 'html', ci=TRUE, ci.level=0.95, single.row = T)

Regression ANOVA

\(R^2 =\) 0.62

kable(anova(reg.Jan.noVax), digits=3, format = "simple", align = 'c', caption = "January Regression ANOVA")

January Regression ANOVA
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
HouseholdIncome.cat	2	61029996.6	30514998.3	87.020	0.000
CityRaceEthnicty.m	3	539142.7	179714.2	0.512	0.674
CityAge.m	2	3066651.9	1533325.9	4.373	0.015
Residuals	113	39625368.5	350667.0	NA	NA

August Incident Rates:

Vaccination Only

Regression Model

reg.Aug.vax <- lm(Aug.IR ~ Aug.Vax.cat)

Regression summary

kable(summary(reg.Aug.vax)$coef, digits=3, format = "simple", align = 'c', caption = "August Regression")

August Regression
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	440.459	21.876	20.134	0.000
Aug.Vax.catHigh	-83.415	38.752	-2.153	0.033
Aug.Vax.catLow	66.055	38.309	1.724	0.087

Regression ANOVA

\(R^2 =\) 0.09

kable(anova(reg.Aug.vax), digits=3, format = "simple", align = 'c', caption = "August Regression ANOVA")

August Regression ANOVA
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Aug.Vax.cat	2	331100.9	165550.46	5.579	0.005
Residuals	118	3501222.4	29671.38	NA	NA

Full Model

Regression Model

reg.Aug <- lm(Aug.IR ~ Aug.Vax.cat + Jan.IR.cat + HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)

Regression summary

kable(summary(reg.Aug)$coef, digits=3, format = "simple", align = 'c', caption = "August Regression")

August Regression
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	303.255	114.684	2.644	0.009
Aug.Vax.catHigh	-17.502	41.118	-0.426	0.671
Aug.Vax.catLow	54.124	40.503	1.336	0.184
Jan.IR.catHigh	-23.978	47.079	-0.509	0.612
Jan.IR.catLow	-232.855	49.476	-4.706	0.000
HouseholdIncome.catHigh	58.249	49.138	1.185	0.238
HouseholdIncome.catLow	-64.624	48.905	-1.321	0.189
CityRaceEthnicty.mCityHispanic.cat	70.693	25.686	2.752	0.007
CityRaceEthnicty.mCityBlack.cat	91.902	44.029	2.087	0.039
CityRaceEthnicty.mCityAsian.cat	-6.602	20.433	-0.323	0.747
CityAge.mCityAgeOld.cat	-72.455	39.920	-1.815	0.072
CityAge.mCityAgeYoung.cat	50.944	46.605	1.093	0.277

Regression ANOVA

\(R^2 =\) 0.3

kable(anova(reg.Aug), digits=3, format = "simple", align = 'c', caption = "August Regression ANOVA")

August Regression ANOVA
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Aug.Vax.cat	2	331100.91	165550.46	6.767	0.002
Jan.IR.cat	2	332743.39	166371.70	6.801	0.002
HouseholdIncome.cat	2	90492.19	45246.10	1.850	0.162
CityRaceEthnicty.m	3	277635.17	92545.05	3.783	0.013
CityAge.m	2	133790.88	66895.44	2.734	0.069
Residuals	109	2666560.76	24463.86	NA	NA

August Vaccination Rates:

Regression Model with August Vaccination

reg.Aug.vax <- lm(Aug.Vax.percent ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)

Regression summary

kable(summary(reg.Aug.vax)$coef, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression")

August Vaccination Regression
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	64.568	7.455	8.661	0.000
Jan.IR.catHigh	-1.109	3.009	-0.369	0.713
Jan.IR.catLow	-0.312	3.090	-0.101	0.920
HouseholdIncome.catHigh	3.086	3.189	0.968	0.335
HouseholdIncome.catLow	0.010	3.169	0.003	0.998
CityRaceEthnicty.mCityHispanic.cat	4.182	1.639	2.551	0.012
CityRaceEthnicty.mCityBlack.cat	1.383	2.857	0.484	0.629
CityRaceEthnicty.mCityAsian.cat	3.774	1.268	2.976	0.004
CityAge.mCityAgeOld.cat	-2.359	2.567	-0.919	0.360
CityAge.mCityAgeYoung.cat	-6.221	2.890	-2.153	0.034

Regression ANOVA

kable(anova(reg.Aug.vax), digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression ANOVA")

August Vaccination Regression ANOVA
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Jan.IR.cat	2	1498.373	749.186	7.238	0.001
HouseholdIncome.cat	2	906.977	453.489	4.381	0.015
CityRaceEthnicty.m	3	1844.707	614.902	5.941	0.001
CityAge.m	2	509.743	254.872	2.462	0.090
Residuals	111	11488.636	103.501	NA	NA

Multnomial Regression Model with August Vaccination

reg.Aug.vax <- lm(Aug.Vax.percent ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
reg.Aug.vax <- multinom(Aug.Vax.cat ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)

## # weights:  33 (20 variable)
## initial  value 132.932087 
## iter  10 value 88.707311
## iter  20 value 88.244184
## final  value 88.244116 
## converged

Regression summary for High Vaccination Rates

coef.high <- summary(reg.Aug.vax)$coefficients["High",]
OR.high <- exp(coef.high)
se.high <- summary(reg.Aug.vax)$standard.errors["High",]
z.high <- coef.high/se.high
p.value.high <- (1 - pnorm(abs(z.high), 0, 1)) * 2
r.high <- data.frame(OR.high, coef.high, se.high, z.high, p.value.high)
kable(r.high, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression for High Vax Rates")

August Vaccination Regression for High Vax Rates
	OR.high	coef.high	se.high	z.high	p.value.high
(Intercept)	0.236	-1.442	2.511	-0.574	0.566
Jan.IR.catHigh	0.094	-2.361	1.365	-1.730	0.084
Jan.IR.catLow	6.029	1.797	0.746	2.409	0.016
HouseholdIncome.catHigh	0.500	-0.694	0.851	-0.815	0.415
HouseholdIncome.catLow	3.031	1.109	1.134	0.978	0.328
CityRaceEthnicty.mCityHispanic.cat	1.596	0.468	0.414	1.128	0.259
CityRaceEthnicty.mCityBlack.cat	0.323	-1.130	1.212	-0.932	0.351
CityRaceEthnicty.mCityAsian.cat	1.048	0.047	0.364	0.128	0.898
CityAge.mCityAgeOld.cat	2.395	0.873	0.655	1.333	0.182
CityAge.mCityAgeYoung.cat	0.659	-0.417	1.266	-0.329	0.742

Regression summary for Low Vaccination Rates

coef.low <- summary(reg.Aug.vax)$coefficients["Low",]
OR.low <- exp(coef.low)
se.low <- summary(reg.Aug.vax)$standard.errors["Low",]
z.low <- coef.low/se.low
p.value.low <- (1 - pnorm(abs(z.low), 0, 1)) * 2
r.low <- data.frame(OR.low, coef.low, se.low, z.low, p.value.low)
kable(r.low, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression for Low Vax Rates")

August Vaccination Regression for Low Vax Rates
	OR.low	coef.low	se.low	z.low	p.value.low
(Intercept)	0.240	-1.427	2.132	-0.669	0.503
Jan.IR.catHigh	0.411	-0.888	0.937	-0.948	0.343
Jan.IR.catLow	4.020	1.391	0.988	1.409	0.159
HouseholdIncome.catHigh	0.517	-0.660	0.971	-0.680	0.496
HouseholdIncome.catLow	0.866	-0.144	0.959	-0.150	0.881
CityRaceEthnicty.mCityHispanic.cat	0.517	-0.660	0.516	-1.280	0.201
CityRaceEthnicty.mCityBlack.cat	1.252	0.225	0.679	0.331	0.741
CityRaceEthnicty.mCityAsian.cat	0.209	-1.566	0.508	-3.080	0.002
CityAge.mCityAgeOld.cat	2.693	0.991	0.768	1.289	0.197
CityAge.mCityAgeYoung.cat	7.712	2.043	0.820	2.491	0.013

Model params over time

Full model fit by date
{r} adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat
No vaccination lag time in this example

Vaccination rates

Modeled as a continuous variable. Plots include dates > 04/01/2021

d <- readRDS("data/analysis_data.rds") 
lac_summarytable_dash <- read.csv("data/LA_County_Covid19_CSA_14day_case_death_table_Clean.csv")

# data clean ----
dout <- d %>% 
  mutate(Date.x = ymd(Date.x), 
         Vax.percent = dose1_all_c_prcent, 
         Jan.IR =  as.numeric(adj_case_14day_rate.y),
         Jan.IR.cat =  relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
                                                ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
                                                       "Med"))), ref="Low")) %>%
  dplyr::select(Date.x, city, adj_case_14day_rate.x, 
                Vax.percent, Jan.IR.cat, 
                HouseholdIncome.cat, CityHispanic.cat, CityBlack.cat, 
                CityAsian.cat, CityAgeOld.cat, CityAgeYoung.cat)




test <- dout %>% 
  # dplyr::filter(Date.x > "2021-03-01") %>% 
  tidyr::nest(data = -Date.x) %>% 
  dplyr::mutate(fit = purrr::map(data, ~lm(adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat, data = .x)),
                tidied = purrr::map(fit, ~broom::tidy(.x), quality = purrr::map(fit, ~glance(.x)))) %>% 
  dplyr::select(-data, -fit) %>% 
  tidyr::unnest(tidied) %>% 
  # dplyr::filter(grepl("Vax", term)) %>%
  arrange(Date.x)



# Vaccination rates ---- 
# after april 2021 (vaccination rates too low prior)
p1 <- ggplot(test %>% filter(grepl("Vax.percent", term), Date.x > "2021-04-01"), aes(Date.x, estimate)) + 
  geom_point() + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
  # geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County", 
                      Date > "2021-04-01") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Race

p1 <- ggplot(test %>% filter(term %in% c("CityAsian.catHigh", "CityBlack.catHigh", "CityHispanic.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Household income

p1 <- ggplot(test %>% filter(term %in% c("HouseholdIncome.catHigh", "HouseholdIncome.catLow")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.2))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

Age category

p1 <- ggplot(test %>% filter(term %in% c("CityAgeOld.catHigh", "CityAgeYoung.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

January incidence rates

p1 <- ggplot(test %>% filter(term %in% c("Jan.IR.catHigh", "Jan.IR.catMed")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
# ggsave(p1, width = 10, height = 8, file = "~/Desktop/test1.png")
p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

Model params over time (vax lag)

Vaccination rates

Modeled as a continuous variable. Plots include dates > 04/01/2021

d <- readRDS("data/analysis_data_14vax_lag.rds") 
lac_summarytable_dash <- read.csv("data/LA_County_Covid19_CSA_14day_case_death_table_Clean.csv")

# data clean ----
dout <- d %>% 
  mutate(Date.x = ymd(Date.x), 
         Vax.percent = dose1_all_c_prcent, 
         Jan.IR =  as.numeric(adj_case_14day_rate.y),
         Jan.IR.cat =  relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
                                                ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
                                                       "Med"))), ref="Low")) %>%
  dplyr::select(Date.x, city, adj_case_14day_rate.x, 
                Vax.percent, Jan.IR.cat, 
                HouseholdIncome.cat, CityHispanic.cat, CityBlack.cat, 
                CityAsian.cat, CityAgeOld.cat, CityAgeYoung.cat)




test <- dout %>% 
  # dplyr::filter(Date.x > "2021-03-01") %>% 
  tidyr::nest(data = -Date.x) %>% 
  dplyr::mutate(fit = purrr::map(data, ~lm(adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat, data = .x)),
                tidied = purrr::map(fit, ~broom::tidy(.x), quality = purrr::map(fit, ~glance(.x)))) %>% 
  dplyr::select(-data, -fit) %>% 
  tidyr::unnest(tidied) %>% 
  # dplyr::filter(grepl("Vax", term)) %>%
  arrange(Date.x)



# Vaccination rates ---- 
# after april 2021 (vaccination rates too low prior)
p1 <- ggplot(test %>% filter(grepl("Vax.percent", term), Date.x > "2021-04-01"), aes(Date.x, estimate)) + 
  geom_point() + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
  # geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County", 
                      Date > "2021-04-01") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Race

p1 <- ggplot(test %>% filter(term %in% c("CityAsian.catHigh", "CityBlack.catHigh", "CityHispanic.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Household income

p1 <- ggplot(test %>% filter(term %in% c("HouseholdIncome.catHigh", "HouseholdIncome.catLow")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.2))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

Age category

p1 <- ggplot(test %>% filter(term %in% c("CityAgeOld.catHigh", "CityAgeYoung.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

January incidence rates

p1 <- ggplot(test %>% filter(term %in% c("Jan.IR.catHigh", "Jan.IR.catMed")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
# ggsave(p1, width = 10, height = 8, file = "~/Desktop/test1.png")
p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

COVID-19 LA County City-level Analysis

David Conti

28 November 2021

The Data.

Hypotheses:

Variable Distributions and Coding

Cities

Incidence Rates

January Peak: 2021-01-07

August Peak: 2021-08-18

Percent Vaccinated in Each City

Race/Ethnicity

Proportion of Hispanic Individuals

Proportion of Black Individuals

Proportion of Asian Individuals

Age

Household Income

Variables for Analysis

Distributions

Correlation

Pairwise Trend Plots

January IR vs. Proportion of Hispanic Individuals

January IR vs. Household Income

August IR vs. January IR

August IR vs. Proportion of Hispanic Individuals

August IR vs. Household Income

August IR vs. August Percent Vaccinated

January Incident Rates:

Regression Model

Regression summary

Regression ANOVA

August Incident Rates:

Vaccination Only

Regression Model

Regression summary

Regression ANOVA

Full Model

Regression Model

Regression summary

Regression ANOVA

August Vaccination Rates:

Regression Model with August Vaccination

Regression summary

Regression ANOVA

Multnomial Regression Model with August Vaccination

Regression summary for High Vaccination Rates

Regression summary for Low Vaccination Rates

Model params over time

Vaccination rates

Race

Household income

Age category

January incidence rates

Model params over time (vax lag)

Vaccination rates

Race

Household income

Age category

January incidence rates